ARIMA

Ejercicio de Simulación Modelos ARIMA

Vamos a ver que existen varios escenarios de raíces unitarias en modelos ARIMA

library(urca)
library(forecast)
library(tseries)
library(lmtest)
library(uroot)
library(fUnitRoots)
library(aTSA)


######Ejercicios de Simulación#####
###################################

##Caminata Aleatoria con y sin drift
set.seed(154) 
w = rnorm(200); x = cumsum(w) 
wd = w +.2; xd = cumsum(wd)
plot.ts(xd, ylim=c(-5,55), main="Caminata Aletoria", ylab='')
lines(x, col=4); abline(h=0, col=4, lty=2); abline(a=0, b=.2, lty=2)

##### Raíz unitaria con componentes ARMA
Tlength=200
a0=3
a1=0.5
tiempo=seq(1:Tlength)
xt=a1*tiempo
tendencia=a0+a1*tiempo
drift=a1*rep(1,Tlength)
arimaej=arima.sim(list(order = c(1,0,1), ar = 0.7,ma=0.6), n = Tlength)

plot(arimaej)

caminata=as.ts(cumsum(arimaej))
drift=as.ts(cumsum(arimaej+a0))


#x11()
#par(mfrow=c(2,1))
plot(caminata)

plot(drift)   ###Caminata Aleatoria Con Drift

linear=as.ts(cumsum(arimaej+a0+xt))

plot(linear)  ###Caminata Aleatoria alrededor de una función cuadrática

auto.arima(arimaej,d=0,D=0,max.p=20,max.q=0,start.p=0, start.q=0,seasonal=FALSE,max.order=20,stationary=TRUE,ic="aic",stepwise=FALSE,allowmean = TRUE)

Series: arimaej 
ARIMA(12,0,0) with zero mean 

Coefficients:
         ar1      ar2     ar3      ar4     ar5      ar6     ar7      ar8     ar9     ar10    ar11     ar12
      1.2429  -0.7827  0.5695  -0.4683  0.4649  -0.3968  0.3822  -0.3601  0.2715  -0.1908  0.3065  -0.1350
s.e.  0.0700   0.1109  0.1238   0.1293  0.1313   0.1329  0.1322   0.1308  0.1295   0.1262  0.1145   0.0719

sigma^2 = 0.9835:  log likelihood = -277.35
AIC=580.71   AICc=582.67   BIC=623.59

ar(arimaej)

Call:
ar(x = arimaej)

Coefficients:
      1        2        3        4        5        6        7        8        9  
 1.2522  -0.7899   0.5314  -0.3955   0.3921  -0.3277   0.3174  -0.2983   0.2001  

Order selected 9  sigma^2 estimated as  1.043

fUnitRoots::adfTest(arimaej,lags = 10,type='nc')

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 10
  STATISTIC:
    Dickey-Fuller: -1.6543
  P VALUE:
    0.09434 

Description:
 Wed Aug 28 12:56:56 2024 by user: 

aTSA::adf.test(arimaej,nlag = 11)

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
      lag   ADF p.value
 [1,]   0 -4.02  0.0100
 [2,]   1 -5.45  0.0100
 [3,]   2 -3.81  0.0100
 [4,]   3 -3.77  0.0100
 [5,]   4 -3.04  0.0100
 [6,]   5 -3.15  0.0100
 [7,]   6 -2.72  0.0100
 [8,]   7 -2.78  0.0100
 [9,]   8 -2.18  0.0298
[10,]   9 -1.96  0.0493
[11,]  10 -1.65  0.0943
Type 2: with drift no trend 
      lag   ADF p.value
 [1,]   0 -4.13  0.0100
 [2,]   1 -5.68  0.0100
 [3,]   2 -3.97  0.0100
 [4,]   3 -3.90  0.0100
 [5,]   4 -3.14  0.0264
 [6,]   5 -3.30  0.0177
 [7,]   6 -2.84  0.0584
 [8,]   7 -2.92  0.0468
 [9,]   8 -2.27  0.2221
[10,]   9 -2.02  0.3185
[11,]  10 -1.70  0.4464
Type 3: with drift and trend 
      lag   ADF p.value
 [1,]   0 -4.55  0.0100
 [2,]   1 -6.26  0.0100
 [3,]   2 -4.50  0.0100
 [4,]   3 -4.35  0.0100
 [5,]   4 -3.58  0.0362
 [6,]   5 -3.87  0.0166
 [7,]   6 -3.37  0.0605
 [8,]   7 -3.52  0.0419
 [9,]   8 -2.81  0.2342
[10,]   9 -2.56  0.3417
[11,]  10 -2.28  0.4566
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

####No hay presencia de Raíz Unitaria
tseries::adf.test(drift,k = 10)   ####Note que esta función no trabaja bien, favor no usar

    Augmented Dickey-Fuller Test

data:  drift
Dickey-Fuller = -1.9002, Lag order = 10, p-value = 0.6179
alternative hypothesis: stationary

fUnitRoots::adfTest(drift,lags = 10,type='c')

Warning in fUnitRoots::adfTest(drift, lags = 10, type = "c") :
  p-value greater than printed p-value

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 10
  STATISTIC:
    Dickey-Fuller: 1.1031
  P VALUE:
    0.99 

Description:
 Wed Aug 28 12:56:57 2024 by user: 

aTSA::adf.test(drift,nlag = 11)

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
      lag   ADF p.value
 [1,]   0 15.74   0.990
 [2,]   1  3.13   0.990
 [3,]   2  4.20   0.990
 [4,]   3  2.89   0.990
 [5,]   4  2.94   0.990
 [6,]   5  2.57   0.990
 [7,]   6  2.64   0.990
 [8,]   7  2.38   0.990
 [9,]   8  2.41   0.990
[10,]   9  1.89   0.985
[11,]  10  1.75   0.979
Type 2: with drift no trend 
      lag  ADF p.value
 [1,]   0 4.76    0.99
 [2,]   1 1.58    0.99
 [3,]   2 2.13    0.99
 [4,]   3 1.45    0.99
 [5,]   4 1.53    0.99
 [6,]   5 1.48    0.99
 [7,]   6 1.52    0.99
 [8,]   7 1.44    0.99
 [9,]   8 1.46    0.99
[10,]   9 1.16    0.99
[11,]  10 1.10    0.99
Type 3: with drift and trend 
      lag    ADF p.value
 [1,]   0 -0.525   0.980
 [2,]   1 -1.375   0.837
 [3,]   2 -1.080   0.923
 [4,]   3 -1.174   0.908
 [5,]   4 -1.139   0.914
 [6,]   5 -1.439   0.811
 [7,]   6 -1.327   0.858
 [8,]   7 -1.510   0.781
 [9,]   8 -1.420   0.819
[10,]   9 -1.699   0.701
[11,]  10 -1.900   0.616
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

###Otro Ejemplo

estacionario=arima.sim(list(order = c(1,0,1), ar = 0.7,ma=0.6), n = Tlength)
trend.estacionario=tendencia+estacionario
#x11()
plot(trend.estacionario)

fUnitRoots::adfTest(trend.estacionario)

Warning in fUnitRoots::adfTest(trend.estacionario) :
  p-value greater than printed p-value

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 1
  STATISTIC:
    Dickey-Fuller: 2.8525
  P VALUE:
    0.99 

Description:
 Wed Aug 28 12:56:57 2024 by user: 

fUnitRoots::adfTest(trend.estacionario,lags=1,type='ct')   ###Cambie 11 rezago hasta 

Warning in fUnitRoots::adfTest(trend.estacionario, lags = 1, type = "ct") :
  p-value smaller than printed p-value

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 1
  STATISTIC:
    Dickey-Fuller: -6.3342
  P VALUE:
    0.01 

Description:
 Wed Aug 28 12:56:57 2024 by user: 

aTSA::adf.test(trend.estacionario,nlag = 2)

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
     lag  ADF p.value
[1,]   0 4.66    0.99
[2,]   1 2.85    0.99
Type 2: with drift no trend 
     lag     ADF p.value
[1,]   0  0.0188   0.957
[2,]   1 -0.2507   0.924
Type 3: with drift and trend 
     lag   ADF p.value
[1,]   0 -4.00  0.0104
[2,]   1 -6.33  0.0100
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

Tlength=200

arimaej_raiz_unit=arima.sim(list(order = c(1,1,1), ar = 0.7,ma=0.6), n = Tlength)
plot(arimaej_raiz_unit)

acf(arimaej_raiz_unit)

pacf(arimaej_raiz_unit)

stats::ar(arimaej_raiz_unit)### Permite Seleccionar el lag

Call:
stats::ar(x = arimaej_raiz_unit)

Coefficients:
      1        2        3  
 1.3077  -0.1686  -0.1512  

Order selected 3  sigma^2 estimated as  7.736

fUnitRoots::adfTest(arimaej_raiz_unit,lags =1 ,type='nc')  

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 1
  STATISTIC:
    Dickey-Fuller: -1.2511
  P VALUE:
    0.2165 

Description:
 Wed Aug 28 12:57:04 2024 by user: 

fUnitRoots::adfTest(arimaej_raiz_unit,type='nc') 

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 1
  STATISTIC:
    Dickey-Fuller: -1.2511
  P VALUE:
    0.2165 

Description:
 Wed Aug 28 12:57:04 2024 by user: 

####No hay presencia de Raíz Unitaria
aTSA::adf.test(arimaej_raiz_unit,nlag=12)   ####Note que hay 

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
      lag    ADF p.value
 [1,]   0 -0.197   0.587
 [2,]   1 -1.251   0.231
 [3,]   2 -0.789   0.396
 [4,]   3 -1.074   0.294
 [5,]   4 -0.841   0.378
 [6,]   5 -0.912   0.352
 [7,]   6 -0.841   0.378
 [8,]   7 -0.975   0.330
 [9,]   8 -0.961   0.335
[10,]   9 -0.799   0.393
[11,]  10 -0.796   0.394
[12,]  11 -0.846   0.376
Type 2: with drift no trend 
      lag   ADF p.value
 [1,]   0 -1.16   0.641
 [2,]   1 -2.11   0.284
 [3,]   2 -1.60   0.484
 [4,]   3 -1.99   0.330
 [5,]   4 -1.63   0.474
 [6,]   5 -1.69   0.450
 [7,]   6 -1.65   0.464
 [8,]   7 -1.75   0.425
 [9,]   8 -1.78   0.416
[10,]   9 -1.71   0.441
[11,]  10 -1.66   0.459
[12,]  11 -1.76   0.423
Type 3: with drift and trend 
      lag    ADF p.value
 [1,]   0  0.203   0.990
 [2,]   1 -2.180   0.499
 [3,]   2 -1.331   0.856
 [4,]   3 -1.908   0.613
 [5,]   4 -1.455   0.804
 [6,]   5 -1.568   0.756
 [7,]   6 -1.469   0.798
 [8,]   7 -1.671   0.713
 [9,]   8 -1.664   0.716
[10,]   9 -1.459   0.802
[11,]  10 -1.447   0.807
[12,]  11 -1.539   0.769
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

######
prueba_df=ur.df(arimaej_raiz_unit,type="none",lags=5)
summary(prueba_df)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.41389 -0.72024  0.00503  0.62760  2.77070 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
z.lag.1     -0.001714   0.001879  -0.912   0.3628    
z.diff.lag1  1.298787   0.072602  17.889  < 2e-16 ***
z.diff.lag2 -0.802310   0.117554  -6.825 1.16e-10 ***
z.diff.lag3  0.533406   0.125282   4.258 3.25e-05 ***
z.diff.lag4 -0.272105   0.117722  -2.311   0.0219 *  
z.diff.lag5  0.065239   0.073395   0.889   0.3752    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.004 on 189 degrees of freedom
Multiple R-squared:  0.7757,    Adjusted R-squared:  0.7686 
F-statistic: 108.9 on 6 and 189 DF,  p-value: < 2.2e-16


Value of test-statistic is: -0.9123 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

Serie de Pasajeros

Vamos a Analizar la Serie de Pasajeros.Iniciamos con las gráficas y transformación de Box-Cox

[1] 4.102259e-05
[1] 0.2
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec
1949 2.548535 2.561426 2.588461 2.582990 2.567558 2.593773 2.615143 2.615143 2.595510 2.563492 2.529884 2.561426
1950 2.555089 2.577352 2.603953 2.593773 2.575434 2.616686 2.646282 2.646282 2.629992 2.590249 2.552929 2.602296
1951 2.610433 2.618215 2.656336 2.636968 2.648852 2.656336 2.680161 2.680161 2.663501 2.635595 2.612017 2.641022
1952 2.647572 2.658759 2.673698 2.659957 2.662328 2.699069 2.709945 2.720109 2.690390 2.671486 2.648852 2.674793
1953 2.676962 2.676962 2.715111 2.714262 2.709067 2.720927 2.737151 2.742897 2.715955 2.692360 2.658759 2.682259
1954 2.685357 2.668111 2.714262 2.707296 2.713408 2.737151 2.762643 2.756996 2.733443 2.709067 2.684331 2.709067
1955 2.720109 2.712549 2.739332 2.740768 2.741481 2.770427 2.796406 2.787934 2.768668 2.744300 2.715955 2.747066
1956 2.751119 2.746379 2.771588 2.769257 2.772164 2.801154 2.818211 2.814887 2.791987 2.765084 2.742191 2.765084
1957 2.770427 2.762027 2.792485 2.788447 2.791987 2.821853 2.837961 2.838662 2.814466 2.787934 2.764478 2.782161
1958 2.784288 2.772164 2.795437 2.788447 2.795922 2.826939 2.846793 2.851301 2.814466 2.793969 2.767483 2.782696
1959 2.794460 2.785340 2.815307 2.811044 2.821052 2.840400 2.864196 2.867286 2.837255 2.815726 2.795437 2.814887
1960 2.819842 2.808860 2.820650 2.836545 2.840400 2.860440 2.883580 2.879651 2.852247 2.836545 2.808419 2.825783

Prueba de Raíz Unitaria

Ahora avanzamos en el setido de verificar si la serie muestra la presencia de una o varias raíces unitarias

acf(logAirP,lag.max = 60)

pacf(logAirP)

ar(logAirP) ##Selecciona un modelo AR usando el crierio de Akaike  a la serie Aaa

Call:
ar(x = logAirP)

Coefficients:
      1        2        3        4        5        6        7        8        9       10       11       12       13  
 1.0013  -0.0959   0.0329  -0.0252   0.0231   0.0005  -0.0157  -0.0762   0.1071  -0.0236   0.0675   0.4536  -0.4854  

Order selected 13  sigma^2 estimated as  0.01336

aTSA::adf.test(logAirP,nlag = 14)

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
      lag    ADF p.value
 [1,]   0  0.913   0.902
 [2,]   1  0.674   0.836
 [3,]   2  0.796   0.871
 [4,]   3  0.936   0.905
 [5,]   4  1.510   0.966
 [6,]   5  1.382   0.957
 [7,]   6  1.451   0.962
 [8,]   7  1.847   0.983
 [9,]   8  3.491   0.990
[10,]   9  4.174   0.990
[11,]  10  8.176   0.990
[12,]  11 11.250   0.990
[13,]  12  3.787   0.990
[14,]  13  2.483   0.990
Type 2: with drift no trend 
      lag    ADF p.value
 [1,]   0 -1.816   0.400
 [2,]   1 -2.018   0.321
 [3,]   2 -1.650   0.465
 [4,]   3 -1.609   0.481
 [5,]   4 -1.288   0.595
 [6,]   5 -1.108   0.659
 [7,]   6 -0.923   0.723
 [8,]   7 -0.807   0.764
 [9,]   8 -0.720   0.795
[10,]   9 -0.859   0.746
[11,]  10 -1.578   0.493
[12,]  11 -2.499   0.134
[13,]  12 -1.952   0.347
[14,]  13 -1.717   0.439
Type 3: with drift and trend 
      lag   ADF p.value
 [1,]   0 -4.85  0.0100
 [2,]   1 -7.00  0.0100
 [3,]   2 -6.71  0.0100
 [4,]   3 -7.13  0.0100
 [5,]   4 -5.66  0.0100
 [6,]   5 -6.42  0.0100
 [7,]   6 -6.88  0.0100
 [8,]   7 -6.28  0.0100
 [9,]   8 -3.62  0.0341
[10,]   9 -2.98  0.1692
[11,]  10 -1.32  0.8577
[12,]  11 -0.44  0.9833
[13,]  12 -1.53  0.7696
[14,]  13 -2.15  0.5109
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

summary(ur.df(logAirP,type="none",lags = 13))

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.090906 -0.022763 -0.002391  0.022181  0.107123 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
z.lag.1       0.003636   0.001464   2.483 0.014448 *  
z.diff.lag1  -0.379495   0.088572  -4.285  3.8e-05 ***
z.diff.lag2  -0.218029   0.071925  -3.031 0.003003 ** 
z.diff.lag3  -0.137513   0.074265  -1.852 0.066618 .  
z.diff.lag4  -0.227017   0.073095  -3.106 0.002386 ** 
z.diff.lag5  -0.098892   0.075307  -1.313 0.191712    
z.diff.lag6  -0.193729   0.071366  -2.715 0.007649 ** 
z.diff.lag7  -0.153409   0.072344  -2.121 0.036090 *  
z.diff.lag8  -0.264234   0.072215  -3.659 0.000382 ***
z.diff.lag9  -0.100257   0.075728  -1.324 0.188138    
z.diff.lag10 -0.204136   0.072692  -2.808 0.005847 ** 
z.diff.lag11 -0.088297   0.073995  -1.193 0.235190    
z.diff.lag12  0.694627   0.072012   9.646  < 2e-16 ***
z.diff.lag13  0.311279   0.088880   3.502 0.000656 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04121 on 116 degrees of freedom
Multiple R-squared:  0.8701,    Adjusted R-squared:  0.8544 
F-statistic:  55.5 on 14 and 116 DF,  p-value: < 2.2e-16


Value of test-statistic is: 2.4833 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

summary(ur.df(logAirP,type="drift",lags = 13))

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.097948 -0.021365 -0.003174  0.022461  0.101422 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.120185   0.056858   2.114 0.036696 *  
z.lag.1      -0.016718   0.009737  -1.717 0.088668 .  
z.diff.lag1  -0.401952   0.087921  -4.572 1.23e-05 ***
z.diff.lag2  -0.253986   0.072886  -3.485 0.000698 ***
z.diff.lag3  -0.180754   0.075984  -2.379 0.019015 *  
z.diff.lag4  -0.266718   0.074435  -3.583 0.000499 ***
z.diff.lag5  -0.143868   0.077196  -1.864 0.064919 .  
z.diff.lag6  -0.234797   0.072957  -3.218 0.001676 ** 
z.diff.lag7  -0.199130   0.074495  -2.673 0.008609 ** 
z.diff.lag8  -0.310371   0.074432  -4.170 5.94e-05 ***
z.diff.lag9  -0.155320   0.079037  -1.965 0.051808 .  
z.diff.lag10 -0.255366   0.075618  -3.377 0.001000 ** 
z.diff.lag11 -0.143773   0.077492  -1.855 0.066112 .  
z.diff.lag12  0.641419   0.075291   8.519 7.14e-14 ***
z.diff.lag13  0.290254   0.088144   3.293 0.001318 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04061 on 115 degrees of freedom
Multiple R-squared:  0.874, Adjusted R-squared:  0.8587 
F-statistic: 56.98 on 14 and 115 DF,  p-value: < 2.2e-16


Value of test-statistic is: -1.717 5.4095 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

######SE debe chequear si hay que diferenciar de nuevo!!!
dlogAirPass=diff(logAirP)
plot(dlogAirPass)

#####Transformación requrida para los datos(transformación Box-Cox y diferencia ordinaria)
pacf(dlogAirPass)

acf(dlogAirPass,lag.max = 48)

ar(dlogAirPass)

Call:
ar(x = dlogAirPass)

Coefficients:
      1        2        3        4        5        6        7        8        9       10       11       12       13  
-0.0993  -0.1226  -0.2824  -0.2875  -0.1184  -0.2406  -0.1829  -0.2971  -0.1040  -0.2627  -0.1259   0.5734   0.0179  
     14       15  
-0.1667   0.1200  

Order selected 15  sigma^2 estimated as  0.002988

aTSA::adf.test(dlogAirPass,nlag = 15)

Augmented Dickey-Fuller Test 
alternative: stationary 
 
Type 1: no drift no trend 
      lag     ADF p.value
 [1,]   0  -9.606  0.0100
 [2,]   1  -8.816  0.0100
 [3,]   2  -7.634  0.0100
 [4,]   3  -8.745  0.0100
 [5,]   4  -6.791  0.0100
 [6,]   5  -6.245  0.0100
 [7,]   6  -6.571  0.0100
 [8,]   7  -9.498  0.0100
 [9,]   8  -8.028  0.0100
[10,]   9 -10.214  0.0100
[11,]  10  -7.327  0.0100
[12,]  11  -1.868  0.0622
[13,]  12  -1.223  0.2402
[14,]  13  -1.283  0.2186
[15,]  14  -0.992  0.3231
Type 2: with drift no trend 
      lag    ADF p.value
 [1,]   0  -9.63  0.0100
 [2,]   1  -8.86  0.0100
 [3,]   2  -7.71  0.0100
 [4,]   3  -8.94  0.0100
 [5,]   4  -6.98  0.0100
 [6,]   5  -6.46  0.0100
 [7,]   6  -6.91  0.0100
 [8,]   7 -10.55  0.0100
 [9,]   8  -9.57  0.0100
[10,]   9 -15.32  0.0100
[11,]  10 -16.05  0.0100
[12,]  11  -4.60  0.0100
[13,]  12  -3.05  0.0352
[14,]  13  -3.22  0.0222
[15,]  14  -2.72  0.0789
Type 3: with drift and trend 
      lag    ADF p.value
 [1,]   0  -9.60  0.0100
 [2,]   1  -8.83  0.0100
 [3,]   2  -7.69  0.0100
 [4,]   3  -8.92  0.0100
 [5,]   4  -6.95  0.0100
 [6,]   5  -6.43  0.0100
 [7,]   6  -6.88  0.0100
 [8,]   7 -10.51  0.0100
 [9,]   8  -9.55  0.0100
[10,]   9 -15.44  0.0100
[11,]  10 -16.57  0.0100
[12,]  11  -4.94  0.0100
[13,]  12  -3.37  0.0628
[14,]  13  -3.56  0.0392
[15,]  14  -3.12  0.1090
---- 
Note: in fact, p.value = 0.01 means p.value <= 0.01 

summary(ur.df(dlogAirPass,type="none",lags = 14))

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.081539 -0.022167  0.001427  0.028987  0.098928 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
z.lag.1      -0.33366    0.33626  -0.992 0.323188    
z.diff.lag1  -0.94573    0.33797  -2.798 0.006042 ** 
z.diff.lag2  -1.06221    0.34615  -3.069 0.002692 ** 
z.diff.lag3  -1.22917    0.34971  -3.515 0.000634 ***
z.diff.lag4  -1.32684    0.32990  -4.022 0.000105 ***
z.diff.lag5  -1.26830    0.30882  -4.107 7.61e-05 ***
z.diff.lag6  -1.34019    0.28801  -4.653 8.95e-06 ***
z.diff.lag7  -1.33665    0.26537  -5.037 1.81e-06 ***
z.diff.lag8  -1.46466    0.24807  -5.904 3.79e-08 ***
z.diff.lag9  -1.41286    0.23488  -6.015 2.26e-08 ***
z.diff.lag10 -1.49206    0.22124  -6.744 6.87e-10 ***
z.diff.lag11 -1.42586    0.20715  -6.883 3.46e-10 ***
z.diff.lag12 -0.60007    0.19564  -3.067 0.002704 ** 
z.diff.lag13 -0.22366    0.15295  -1.462 0.146431    
z.diff.lag14 -0.20779    0.09268  -2.242 0.026919 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04183 on 113 degrees of freedom
Multiple R-squared:  0.9165,    Adjusted R-squared:  0.9054 
F-statistic: 82.71 on 15 and 113 DF,  p-value: < 2.2e-16


Value of test-statistic is: -0.9923 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

monthplot(dlogAirPass)

spectrum(dlogAirPass)

Estimated spectral density
    series: dlogAirPass
    method: Raw Periodogram
    nseasons: 12
    frequency range: (-6,6]

    sort method for the table: decreasing magnitudes

             freq       spec    % Total    Cum. %
    [1,] 2.000000 1.73274552 0.28973097 0.2897310
    [2,] 1.000000 1.61152998 0.26946262 0.5591936
    [3,] 4.000000 0.92289692 0.15431684 0.7135104
    [4,] 5.000000 0.48378566 0.08089341 0.7944038
    [5,] 3.000000 0.46145560 0.07715962 0.8715635
    [6,] 4.166667 0.08107252 0.01355607 0.8851195

Identificación del Modelo Usando ACF y PACF

####Identificación de los Órdenes Autoregresivos
acf(dlogAirPass,lag.max = 60,ci.type='ma')  ###Se requiere un MA de orden muy grande

pacf(dlogAirPass,lag.max = 60)  ####Puede ser un autoregresivo de orden 12

###Arima Automático
modelo.automatico1=auto.arima(dlogAirPass,d=0,D=0,max.p=12,max.q=48,start.p=0, start.q=0,seasonal=FALSE,max.order=12,stationary=TRUE,ic="aicc",stepwise=FALSE,allowmean = TRUE)

Ajuste del Modelo

#####Ajuste del Modelo
####Note que entramos la serie original
library(TSA)

Attaching package: ‘TSA’

The following objects are masked from ‘package:PerformanceAnalytics’:

    kurtosis, skewness

The following object is masked from ‘package:readr’:

    spec

The following object is masked from ‘package:sarima’:

    periodogram

The following objects are masked from ‘package:stats’:

    acf, arima

The following object is masked from ‘package:utils’:

    tar

AjusteArima93=forecast::Arima(AirPassengers,order = c(9,1,3),lambda = 0,include.drift = TRUE)
summary(AjusteArima93)

Series: AirPassengers 
ARIMA(9,1,3) with drift 
Box Cox transformation: lambda= 0 

Coefficients:
          ar1      ar2     ar3      ar4     ar5      ar6      ar7      ar8      ar9     ma1      ma2      ma3
      -0.4424  -0.2037  0.0860  -0.6827  -0.271  -0.1659  -0.1491  -0.6805  -0.2141  0.4821  -0.1477  -0.7510
s.e.   0.1105   0.0758  0.0759   0.0713   0.100   0.0711   0.0750   0.0715   0.1037  0.0821   0.0677   0.0689
       drift
      0.0101
s.e.  0.0009

sigma^2 = 0.004642:  log likelihood = 183.78
AIC=-339.55   AICc=-336.27   BIC=-298.07

Training set error measures:
                     ME     RMSE      MAE        MPE     MAPE     MASE     ACF1
Training set -0.2192766 20.40352 15.27227 -0.2117784 5.395308 0.476807 0.012135

coeftest(AjusteArima93)

z test of coefficients:

         Estimate  Std. Error  z value  Pr(>|z|)    
ar1   -0.44238067  0.11049430  -4.0037 6.237e-05 ***
ar2   -0.20368876  0.07583340  -2.6860  0.007231 ** 
ar3    0.08599219  0.07591330   1.1328  0.257311    
ar4   -0.68273511  0.07128722  -9.5772 < 2.2e-16 ***
ar5   -0.27099739  0.10000994  -2.7097  0.006734 ** 
ar6   -0.16591360  0.07107162  -2.3345  0.019572 *  
ar7   -0.14906334  0.07498938  -1.9878  0.046835 *  
ar8   -0.68051426  0.07152435  -9.5144 < 2.2e-16 ***
ar9   -0.21410059  0.10371315  -2.0644  0.038984 *  
ma1    0.48207457  0.08205276   5.8752 4.224e-09 ***
ma2   -0.14767013  0.06771628  -2.1807  0.029204 *  
ma3   -0.75098362  0.06886761 -10.9047 < 2.2e-16 ***
drift  0.01009809  0.00090294  11.1836 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#####Refinando el modelo
AjusteArima12=forecast::Arima(AirPassengers,order = c(12,1,0),lambda = 0,include.drift = TRUE)
#,fixed=c(NA,NA,0,NA,NA,NA,NA,NA,NA,NA,0,NA,NA)
summary(AjusteArima12)

Series: AirPassengers 
ARIMA(12,1,0) with drift 
Box Cox transformation: lambda= 0 

Coefficients:
          ar1      ar2      ar3      ar4      ar5      ar6     ar7      ar8      ar9     ar10    ar11    ar12
      -0.2212  -0.2835  -0.2550  -0.2973  -0.2299  -0.2622  -0.257  -0.3274  -0.2194  -0.2934  -0.198  0.6020
s.e.   0.0670   0.0679   0.0673   0.0699   0.0679   0.0657   0.066   0.0679   0.0684   0.0681   0.068  0.0683
       drift
      0.0098
s.e.  0.0011

sigma^2 = 0.001807:  log likelihood = 246.66
AIC=-465.31   AICc=-462.03   BIC=-423.83

Training set error measures:
                     ME     RMSE      MAE          MPE     MAPE      MASE       ACF1
Training set -0.1137697 10.82763 8.086473 -0.007108351 3.093656 0.2524632 -0.2063749

coeftest(AjusteArima12)

z test of coefficients:

        Estimate Std. Error z value  Pr(>|z|)    
ar1   -0.2212162  0.0669842 -3.3025 0.0009582 ***
ar2   -0.2834744  0.0679448 -4.1721 3.018e-05 ***
ar3   -0.2549912  0.0672627 -3.7910 0.0001501 ***
ar4   -0.2973403  0.0698594 -4.2563 2.079e-05 ***
ar5   -0.2298510  0.0679217 -3.3841 0.0007142 ***
ar6   -0.2621654  0.0657330 -3.9883 6.654e-05 ***
ar7   -0.2570036  0.0660119 -3.8933 9.889e-05 ***
ar8   -0.3274286  0.0679449 -4.8190 1.443e-06 ***
ar9   -0.2194360  0.0684078 -3.2078 0.0013377 ** 
ar10  -0.2934215  0.0680902 -4.3093 1.638e-05 ***
ar11  -0.1979913  0.0680155 -2.9110 0.0036031 ** 
ar12   0.6020203  0.0682917  8.8154 < 2.2e-16 ***
drift  0.0097712  0.0010737  9.1003 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Análisis de Residuales

residuales=AjusteArima12$residuals
plot(residuales)

#
#plot(SDresiduales)
acf(residuales,lag.max = 48)

pacf(residuales)

#acf(SDresiduales)
#pacf(SDresiduales)


#Test de normalidad
jarque.bera.test(residuales)

    Jarque Bera Test

data:  residuales
X-squared = 7.0597, df = 2, p-value = 0.02931

#Test de autocorrelación
length(residuales)/4

[1] 36

sqrt(length(residuales))

[1] 12

Box.test(residuales, lag =36 , type = "Ljung-Box", fitdf = 12)

    Box-Ljung test

data:  residuales
X-squared = 64.81, df = 24, p-value = 1.298e-05

monthplot(residuales)

spectrum(residuales,spans = c(3,5))

Estimated spectral density
    series: residuales
    method: Smoothed Periodogram
    nseasons: 12
    frequency range: (-6,6]

    sort method for the table: decreasing magnitudes

             freq      spec    % Total     Cum. %
    [1,] 5.583333 0.2454543 0.04105811 0.04105811
    [2,] 5.666667 0.2217983 0.03710108 0.07815919
    [3,] 5.500000 0.2216797 0.03708124 0.11524043
    [4,] 5.416667 0.1780115 0.02977668 0.14501711
    [5,] 5.750000 0.1645269 0.02752105 0.17253816
    [6,] 2.666667 0.1408333 0.02355775 0.19609590

###Estad?sticas CUSUM
res=residuales
cum=cumsum(res)/sd(res)
N=length(res)
cumq=cumsum(res^2)/sum(res^2)
Af=0.948 ###Cuantil del 95% para la estad?stica cusum
co=0.14013####Valor del cuantil aproximado para cusumsq para n/2
####Para el caso de la serie de pasajeros es aprox (144-12)/2=66
LS=Af*sqrt(N)+2*Af*c(1:length(res))/sqrt(N)
LI=-LS
LQS=co+(1:length(res))/N
LQI=-co+(1:length(res))/N
plot(cum,type="l",ylim=c(min(LI),max(LS)),xlab="t",ylab="",main="CUSUM")
lines(LS,type="S",col="red")

lines(LI,type="S",col="red")

#CUSUMSQ
plot(cumq,type="l",xlab="t",ylab="",main="CUSUMSQ")                      
lines(LQS,type="S",col="red")                                                                           

lines(LQI,type="S",col="red")

#####Fase de Pronósticos
pronosticos12=forecast::forecast(AjusteArima12,h=12,level=0.95)
plot(pronosticos12)